classdef h2params < handle
  properties (Constant)
    A=sqrt(3)*pi/18 - log(3)/2;
    %A = 0;
  end
  properties (GetAccess='public', SetAccess='public')
    % time steps
    Tbegin = 0
    Tfinal
    Tsteps

    % background T offset
    TempOffset   = 1

    % number of samples in each dimension
    Nx        = 32
    Ny        = 32
    Nz        = 32

    % physical size in each dimension
    Lx        = 100
    Ly        = 100
    Lz        = 1
    Lz0       = 0

    % physical cell values 
    Cx
    Cy
    Cz

    % physical distance between adj cells
    Dx
    Dy
    Dz

    Dt
  end

  properties
    % control flags
    metricOn = false
    metricOrder = 2

    metricMark = false

    metricEvaluationPoints = [0.025];
    metricLinearFnInterpolation = false;
    
    plotHydro = false
    plotMetric = false

    colormapOffset = 512

    hydroOrder = 0

    hydroRandGradOn   = false;
    hydroRandGradFreq = 0;

    hydroPureWave = false;
    hydroPureWaveFreq = 1;
    hydroPureWaveAmp = 0.01;

    hydroPeak = false;
    hydroPeakAmp = 0.1;
    hydroPeakRelWidth = 0.1;

    hydroVortices = false
    hydroVorticesParams = [] % [x, y, width, strength]

    hydroSaveOutput = false;

    timeOrder = 2;

  end
  methods 
    function obj=h2params(tfinal, TIME_STEPS, X_GRID_SIZE, L, Z_GRID_SIZE, TEMP_OFFSET)
      obj.Tfinal = tfinal;
      obj.Tsteps = TIME_STEPS;

      obj.TempOffset = TEMP_OFFSET;
      
      % number of samples
      obj.Nx = X_GRID_SIZE; 
      obj.Ny = X_GRID_SIZE;
      obj.Nz = Z_GRID_SIZE;


      obj.Lx = L;
      obj.Ly = L;
      obj.Lz = 1/obj.TempOffset; % apparent horizon

      % min Lz:
      ref0 = 0.1;  % reference boundary at T = 1
      % changes in T will move the position of the horizon.
      % this affects integration of tensors.
      % set a reference: when T = 1, integrate from rh = 1 to ref0 \approx 0
      % at a differentt temperature T, choose the bounday integration limit 
      % so as to get the same volume integral. 
      % this is the answer:
      obj.Lz0 = ref0.*(1+(-1).*ref0.^3+TEMP_OFFSET.^3.*ref0.^3).^(-1/3);
      obj.Lz0 = 0;

      % cell values
      obj.Cx = (0:obj.Nx-1)/obj.Nx*obj.Lx;
      obj.Cy = (0:obj.Ny-1)/obj.Ny*obj.Ly;
      %obj.Cz = (0:obj.Nz-1)/obj.Nz*obj.Lz;
      %obj.Cz = chebpoints(obj.Nz, 0, obj.Lz);
      obj.Cz = chebpoints(obj.Nz, obj.Lz0, obj.Lz);

      obj.Dx = obj.Lx / obj.Nx;
      obj.Dy = obj.Ly / obj.Ny;
      obj.Dz = obj.Lz / obj.Nz;

      obj.Dt = obj.Tfinal ./ obj.Tsteps;
    end

  end
end
